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We present a new impurity solver for dynamical mean-field theory based on imaginary-time evo¬ 
lution of matrix product states. This converges the self-consistency loop on the imaginary-frequency 
axis and obtains real-frequency information in a final real-time evolution. Relative to computations 
on the real-frequency axis, required bath sizes are much smaller and less entanglement is generated, 
so much larger systems can be studied. The power of the method is demonstrated by solutions of a 
three band model in the single and two-site dynamical mean-field approximation. Technical issues 
are discussed, including details of the method, efficiency as compared to other matrix product state 
based impurity solvers, bath construction and its relation to real-frequency computations and the 
analytic continuation problem of quantum Monte Carlo, the choice of basis in dynamical cluster 
approximation, and perspectives for off-diagonal hybridization functions. 


I. INTRODUCTION 

Dynamical mean-field theory (DMFT) in its single¬ 
site [1-3] and cluster [4, 5] variants is among the most 
widely employed computational techniques for solving 
quantum many-body problems. At the core of a numer¬ 
ical solution of DMFT is an impurity solver : an algo¬ 
rithm for solving a quantum impurity problem. The most 
prominent examples of impurity solvers are continuous 
time quantum Monte Carlo (CTQMC) [6-8], exact cliag- 
onalization (ED) [9—11], the numerical renormalization 
group (NRG) [12], and the density matrix renormaliza¬ 
tion group (DMRG) [13]. While all methods have their 
strengths, key limitations mean that fundamentally im¬ 
portant classes of problems have not yet been adequately 
addressed. Other recent suggestions for impurity solvers 
[14-19] including in particular the computationally inex¬ 
pensive density matrix embedding theory (DMET) [20], 
are promising but have not been tested in detail. 

CTQMC is widely employed but its application to sit¬ 
uations involving low point symmetry, non-Hubbard in¬ 
teractions or multiple relevant orbitals is limited by the 
fermionic sign problem. Reaching low temperatures be¬ 
comes highly computationally expensive while calculat¬ 
ing real-frequency information requires analytical contin¬ 
uation, a numerically ill-posed procedure fraught with 
practical difficulties. 

ED makes no assumption on the interaction and does 
not have a sign problem. It is limited by the size of the 
Hilbert space that can be studied, meaning in practice 
that it is restricted to a small number of correlated sites 
to which only a small number of bath sites can be at¬ 
tached. Recently, improvements have been achieved by 
considering only restricted subspaces of the Hilbert space 
[21-24], but the size of problem remains a significant lim¬ 
itation. 

NRG converges the DMFT loop on the real-frequency 
axis and very effectively obtains real-frequency informa¬ 
tion in the low frequency limit. Current applications 


have been to relatively small problems (the most recent 
achievement is a solution of the single-site DMFT ap¬ 
proximation to a three-band model [25]) and it remains 
to be seen how far the method can be extended. 

DMRG [26] is a set of algorithms operating on the 
space of matrix product states (MPS) [27]. It has 
been found to be extremely powerful for the calculation 
of ground states of one-dimensional quantum systems 
27, 28]; it was very successfully extended to the calcu¬ 
lation of spectral functions which, in contrast to NRG, 
it obtains with equal resolution across the spectrum, see 
e.g. [29, 30]. In pioneering work the method was applied 
as a DMFT solver by Garcia, Hallberg, and Rozenberg 
[13] and Nishimoto, Gebhard, and Jeckelmann [31] with 
important further work done by these and other authors 
[32-38]. However, the method has not been widely ac¬ 
cepted, perhaps because high-quality data were presented 
only for the single-site approximation to the single-band 
Hubbard model. Recently the method was used to reli¬ 
ably solve a two-site DCA [39], and insights into the en¬ 
tanglement of the impurity problem made it even more 
powerful [40]. In view of these advances, DMRG now is 
a candidate for a highly flexible low-cost impurity solver, 
which can in addition be efficiently employed in the non¬ 
equilibrium formulation of DMFT [40-42] . 

This paper reformulates the DMRG method as an 
impurity solver for DMFT on the imaginary-time axis 
(Sec. II). As we will show, this strongly reduces en¬ 
tanglement and necessary bath sizes. The price to be 
paid is a reduced resolution on the real-frequency axis, 
which we study in detail by comparing with calculations 
that converge the DMFT loop on the real-frequency axis 
(Sec. III). The reformulation allows a much larger class 
of problems to be addressed, including some that are un¬ 
reachable by other methods, due e.g. to the sign problem, 
the size of the correlated cluster or the number of bands. 
We illustrate this with calculations for three band mod¬ 
els in the single-site and two-site DMFT approximation 
(Sec. IV). We append discussions of the optimization of 
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typical DMFT Hamiltonians (Appendix A) and of the 
entanglement in different representations of the DCA in¬ 
cluding a discussion of off-diagonal hybridization func¬ 
tions (Appendix B). 


II. METHOD 

A. Overview: Green’s functions in DMRG 

The computational key challenge in DMFT is the com¬ 
putation of the full frequency dependence of the Green’s 
function of a quantum impurity model involving an essen¬ 
tially arbitrary bath. The “size” (number of correlated 
sites L c ) of the impurity model should be as large as 
possible and the kinds of interaction that can be treated 
should be as general as possible. The Green’s function is 
used in a self-consistency loop, which may require many 
iterations for convergence. The solution should be as in¬ 
expensive as feasible, and must run automatically, with¬ 
out need for manual optimization of parameters or pro¬ 
cedures. In this subsection we present a qualitative dis¬ 
cussion of the issues involved in computing the Green’s 
function using DMRG methods, to motivate the work 
described in detail below. 

Within DMRG one computes Green’s functions by first 
representing the system ground state \E 0 ) as an MPS. 
One then generates a one electron (one hole) excitation 
| i/>q) = d^\E 0 ) (IV’o') = d\E 0 }) by applying a creation (an¬ 
nihilation) operator d ' (d) to \E 0 ). While the state |^) 
is at most as entangled as the ground state \E 0 ) [30], in 
order to compute a Green’s function, one has to perform 
further operations on | ifg). These operations typically 
increase entanglement and by that the bond dimension 
of an MPS, which ultimately limits all computations. 

Let us be more concrete and consider a general MPS 
of bond dimension m for a system with L sites and open 
boundary conditions. Defining A a ', B rTi £ C mxm for i ^ 
1, L and' A ai £ C lxm , B aL £ C mxl , where cq £ {0,t4 
,tl} labels a local basis state of the Hilbert space, any 
MPS can be represented as [27] 

|^mps)= 51 A"E..A°'SB°'+\..B^\ ai ,...,a L ), (1) 

ai,...,a L 

where S = diag(si,..., s m ) is a diagonal matrix, and A ai 
are left-normalized and B ai are right-normalized , respec¬ 
tively: 

A" 1 * A* = I, 55 B^B^ = I. (2) 

CTj 

Here, / are identity matrices. Left- and right- 
normalization make Eq. (1) the Schmidt decomposition 
of It/’MPs) that is associated with partitioning the sys¬ 
tem at bond (1,1 + 1). The bond entanglement entropy 
for the associated reduced density matrix can therefore 


simply be read of from Eq. (1) [27] 

m 

( 3 ) 

I/=l 

When subsequently we refer to an entanglement growth 
associated with repeated operations on |?/ , mps) this im¬ 
plies the need to adjust the bond dimension m such that 
| V’mps) still faithfully represents a physical state. If en¬ 
tanglement in the physical state becomes too large, we 
have to choose m so large that computations with MPS 
become impractical. 

Since the first suggestion for computing spectral func¬ 
tions within DMRG [43], the field has evolved by the 
important development of the correction vector method 
[44, 45] . The subsequent understanding of the connection 
between DMRG and MPS [27] opened the door to many 
further approaches to computing spectral and Green’s 
functions, in particular, time evolution and subsequent 
Fourier transform [46, 47], an improved Lanczos algo¬ 
rithm [ S], and the Chebyshev recursion [29, 30, 49]. All 
of these are formulated for the calculation of spectral 
functions at T = 0, as considered in the present pa¬ 
per, and came at much cheaper computational cost than 
the correction vector method [29, 30]. We note that for 
T > 0, there are perspectives for even more powerful al¬ 
gorithms: it was recently demonstrated that the numer¬ 
ically exact spectral function of a molecule consisting of 
several hundreds of interacting spins could be computed 

[50]. 

These developments (see Appendix C for more details) 
make MPS-based solvers an attractive possibility for dy¬ 
namical mean-field theory. However, the growth of en¬ 
tanglement arising in all calculations of Green’s function 
has limited the systems sizes that have been addressed to 
date. Also, in MPS computations manual adjustments, 
for example choosing optimal broadening [29] or combin¬ 
ing results of different systems sizes [48] , are still common 
practice. In the rest of this section, we show that these 
problems can to a large degree be circumvented by com¬ 
puting Matsubara Green’s functions using imaginary- 
time evolution. The imaginary-time framework naturally 
extends existing approaches based on real-time evolution 
[38, 40], which have recently been shown to provide cur¬ 
rently the most efficient algorithmic approach to compute 
real-frequency spectral functions [30]. 

B. Imaginary-time computation 

The central objects of technical interest in this paper 
are the greater and the lesser correlation functions G^, 
which we define for imaginary time r 

G^(t) = (i’f\e T(H ~ Eo)T \^f), (4a) 

= (i/>$\e^ H ~ E ^}). (4b) 

In the second line, we evaluate G^(r) | r=it and by that 
obtain a correlation function for real time t , which will 
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FIG. 1. (Color online) Panel (a): Imaginary time correlation 
functions G^(t) defined in Eq. (4a) for an impurity model 
arising in the context of the two-site dynamical cluster ap¬ 
proximation to the single-band Hubbard model on the square 
lattice with next-nearest neighbor hopping t'/t = 0.3, half¬ 
band width D = 4t, interaction U = 2.5D and band filling 
n = 0.96 in the paramagnetic phase. See Ref. 51 for defi¬ 
nition of the model and the meaning of the orbital (patch) 
quantum number U K = ±”. The dashed line is obtained us¬ 
ing linear prediction for times tD > 200. Panel (b): Maximal 
bond dimension m of time evolved states. The MPS computa¬ 
tion used a Hamiltonian representation of the impurity model 
with L c = 2 correlated sites and Lb = 14 bath sites. The hy¬ 
bridization function of the impurity model T dlscr was fitted 
using Eq. (12) for /? e ff = 315/D and a = 0. For the ground 
state optimization, we enforced a maximal bond dimension 
of m = 300. The Krylov time evolution used a time step 
of At = 0.1/D and allowed for a maximal global truncation 
error of 10” 4 at each time step adjusting bond dimensions au¬ 
tomatically. This leads to an immediate decay of m at r ~ 0 
from m = 300 down to m ~ 110, as seen in panel (b). We 
used the global SU(2) symmetry of the Hamiltonian to reach 
these low values of the bond dimension. 


be useful later on. The functions carry spin and or¬ 
bital indices associated with the spin and orbital indices 
of the single-particle (hole) excitation | ipo), but these 
indices are not explicitly written here. We will discuss 
the relationship of G^ to the physical Green’s functions 
(which we denote by G) below. 

While it is not essential in principle, we evaluate 
Eq. (4) using a Krylov algorithm [52], which represents 
the time evolution operator in a local Krylov space and is 
able to treat Hamiltonians with long-ranged interactions. 
Before performing a time evolution computation, one has 
to compute the initial state |'i/ , o') using an MPS optimiza¬ 
tion of the ground state. As impurity models come with 


open boundary conditions, this is well suited for DMRG. 
We discuss this optimization for typical DMFT Hamilto¬ 
nians in Appendix Al. 

Figure 1 presents representative results based on pa¬ 
rameters obtained from a two-site DMFT solution of 
the Hubbard model. Figure 1(a) shows the time evo¬ 
lution of G^(t) out to times as long as 350 times the 
basic timescale (inverse half-bandwidth D) of the model, 
which suffices to converge G^(r) to a precision of 5TO” 4 . 
Figure 1(b) demonstrates the key advantage that makes 
this computation possible: the lack of growth of max¬ 
imal bond dimensions m with time of the associated 
imaginary-time evolved states \ip^(r)) = / /,<)_ 

The imaginary time-evolution operator does not cre¬ 
ate entanglement as it projects on the lowly entangled 
ground state. 

Figure 1(a) reveals additional information about the 
nature and rate of convergence of G^ (r). In the insulat¬ 
ing phase, H has a gap and G^ (r) decays exponentially 
irrespective of whether one considers a finite system or 
the thermodynamic limit. In the metallic phase, G^(r) 
decays algebraically in the thermodynamic limit. For a 
finite system though, there always remains a small gap, 
and even though the decay resembles an algebraic decay 
for short times, it always becomes exponential at long 
times. The exponential decay can be exploited to speed 
up computations considerably by a simple technique 
known as linear prediction [30, 47, 53, 54]. This technique 
is an efficient formulation of the fitting problem for the 
ansatz function /(r) = a n e /3riT , a ra ,/3 n £ C,r £ t, 
which can then be used to reliably extrapolate functions 
with an exponentially decaying envelope. This is illus¬ 
trated by the dashed black line in Fig. 1(a), which has 
been fitted to match G^(r) for tD €E [150,200] and was 
then extrapolated for higher times. The solid green line, 
by contrast, is the result of the MPS computation. Agree¬ 
ment can be seen to be perfect. 


C. Physical Green’s functions 

Of particular interest in the rest of this paper are the 
imaginary time Green’s functions G mat (r) defined via 

G<mat ( T ) = —0 (t)G>(t) + 0(-t)G<(t), (5) 

whose Fourier transform gives the Matsubara Green’s 
function 

/ OO 

dr e iw " T G mat (r), (6) 

-OO 

at zero temperature, where u> n = (2 n + l)7r//3 and /3 — > 
oo. We shall also be interested in the retarded real-time 
Green’s function 

G ret (t) = -id(t)(G>(it) + G< (it)), 


(7) 
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from which the retarded frequency-dependent Green’s 
function is obtained as 

/ OO 

dte i{u+l0+)t G iet (t). (8) 

-OO 

In numerical practice, we evaluate the Fourier trans¬ 
forms leading to (6) and (8) approximately as 

/'1’max fO 

G mat («w n ) = — / drG > (r)e“" i + 1 dr G < (r)e i “" t , 

J 0 j— Tmax 

/* ^max 

G ret (w) =-i dt [ G>(it ) + G<(it)]e^, (9) 

Jo 

with cutoff times r max and f max . This approximation 
is only controlled if we are able to reach long enough 
times r max and t max , such that G^(r) and G^(it) have 
converged to zero to any desired accuracy. 

In contrast to a computation on the imaginary axis, 
reaching arbitrarily long times f max on the real axis 
is prohibited by a logarithmic growth of entanglement, 
which comes with a power-law growth of bond dimen¬ 
sions. In addition, finite-size effects are a severe source 
of errors because the long-time behavior is determined 
by the bath size. For a numerically exact computation, 
one has to choose the system large enough to observe ex¬ 
ponential “pseudo-convergence” of G^ (it) to zero before 
finite-size effects are resolved [30]. In the context of the 
present paper, we deal with small system sizes and will 
never observe pseudo-convergence. In particular there is 
no exponential pseudo-convergence, so that linear predic¬ 
tion cannot be employed [30]. Therefore, when comput¬ 
ing the real-frequency spectral function after converging 
the DMFT loop, one has to use the further approxima¬ 
tion of damping the finite-size effects that emerge at long 
times by computing, instead of G ret (w) in (9), 

pOO 2 2 

G T v et (to) = M dt[G > (U) + G < (it)]e iu,t e~ !L ^ i (10) 

Jo 

which yields the broadened spectral function A v (to) = 

-£lm G r v et (ia) = J duo'A(to')e . Instead of a 

Gaussian damping and broadening, one could also use an 
exponential damping leading to Lorentzian broadening, 
which damps out the original time-evolution information 
more strongly, though. 

Before presenting detailed benchmark results for the 
solution of DMFT using imaginary-time evolution of 
MPS, let us clarify which price we have to pay for profit¬ 
ing from the great advantage of not facing entanglement 
growth. We do this by comparing the imaginary-time 
approach (itMPS) to approaches that solve the DMFT 
loop on the real axis. 


III. COMPARISON OF IMAGINARY-AXIS 
WITH REAL-AXIS COMPUTATIONS 

The self-consistency equation in DMFT relates an im¬ 
purity model specified by a hybridization function and a 



FIG. 2. (Color online) Fit of the hybridization function in the 
two-site DCA problem studied in Figs. 1 and 3, but here for 
the case U = 0. The minimization (12) was done using a = 0 
and a frequency grid defined by /3 e ff = 100/D and a cutoff 
frequency of w c = 6 D. Evidently, the quality of the fit does 
not improve any more for L},/L c > 9. 


self energy to a lattice model specified by a lattice Hamil¬ 
tonian and the same self energy. We discuss the issues 
using the example of the dynamical cluster approxima¬ 
tion to the single-band Hubbard model 


GT(z) 


— v_-_ (ii) 

N k ‘§p K z + M - £fc - £k(z) 

[z + p - £k - E K (z) - A K (z)\ 1 =G™ p (z). 


Here E k denotes the single-particle dispersion of the lat¬ 
tice and p is the chemical potential. In the dynamical 
cluster approximation, the Brillouin zone, consisting in 
N momentum vectors k, is covered by N c (for single¬ 
band L c = N c ) equal-area tiles (patches), labelled here 
by Pk and the self energy is taken to be piecewise con¬ 
stant, with being a potentially different function of fre¬ 
quency in each tile. The impurity model is specified 
by the on-site energy ek and the hybridization function 
Ak(z), which is to be determined using the fixed point 
iteration referred to as the DMFT loop. The loop works 
as follows: make an initial guess for Ak(z); compute 
Ek(z) = to + p — ek — A k (z) — [G™ p (z)] _1 , then update 
A k via A k (z) = z + p - ek ~ E k (z) - [G^^z)] -1 and 
repeat this procedure until convergence. 

We discuss two aspects of the comparison of real- 
and imaginary-frequency solutions of the DMFT self- 
consistency equation (11). The first has to do with the 
number of bath sites needed to obtain a solution of the 
self-consistency equation. The second is the accuracy to 
which the spectral functions of physical interest can be 
reproduced. 

The DMFT self-consistency equation (11), defines the 
hybridization function Ak as a continuous function in 
terms of the difference between the computed self en¬ 
ergy and the inverse of the lattice Green’s function. In 
DMRG-type methods, the hybridization function Ak is 
approximated as the hybridization function A^i scr of a 
discrete impurity model, which is the sum of a number 
of poles. If the number of poles is too small one can- 
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not construct a meaningful approximation on the real 
axis [55] and a DMFT loop cannot be converged. For 
this reason DMRG-based solutions of DMFT up to now 
[13, 31-40], all of which were real axis computations, have 
been performed using numbers of bath sites of at least 
Lf) / Lq > 30, and in the case of the single-band Hub¬ 
bard model, even much more L;,/L c g) 120. Use of such 
a large number of bath sites means that with modest 
broadening the hybridization function can be reasonably 
approximated as a continuum, enabling a stable solution 
of Eq. (11). 

By contrast, formulating the problem on the imaginary 
axis (as is typically done in standard ED solvers where 
the number of bath sites is strictly limited) automatically 
smoothens the hybridization function H^l scr and permits 
a stable solution. From the imaginary axis solution one 
must then determine the discrete set of bath parameters 
to represent Ajf CI . This is typically done [9, 11, 56] by 
numerical minimization of a cost function defined as 

, JV £it 

x 2 = Jj- Y.“n a \A K (iu n ) -A^(iu n )\ 2 . (12) 

n— 1 

Here, a defines a weighting function u>~ a . Choosing 
a > 0, e.g. a = 1, attributes more weight to smaller 
frequencies [11, 56, 57], which we find helpful when us¬ 
ing small bath sizes Lb/L c < 5. To define the frequency 
grid for the fit oj n = (2 n + l)n//3 e g, one defines a fic¬ 
titious inverse temperature /3 e ff, which has no physical 
significance. We further employ a cutoff frequency u ) c , 
which implies a finite number Ng t of fitted Matsubara 
frequencies. 

If one tries to define an analogous cost function for the 
real axis, the result is useless as then Ajf CI (uj + i0 + ) 
is a sum of poles, whereas the hybridization function 
Ak{uj + *0+), as encountered in (11), is continuous [55]. 
One can overcome this problem only when using a Lind- 
bladt formalism [58], which increases the complexity of 
the problem substantially. 

The minimization of (12) is done using using standard 
numerical optimization. The optimization in the initial 
DMFT iteration should be done using a global optimiza¬ 
tion scheme [ ], and in subsequent iterations using a lo¬ 

cal optimization scheme (e.g. conjugate gradient), which 
takes as an initial guess for the new bath parameters 
the values of the previous iteration. Figure 2 shows the 
convergence of the fit of the hybridization function with 
the number of bath sites Lt/L c . For Lb/L c = 7 one al¬ 
ready obtains errors as little as ~ 10 -3 and for values 

/L c > 9 the quality of the fit already stops improving. 
It is at this point, where we (and all ED-like techniques) 
face the problem of “analytic continuation” encountered 
in CTQMC, namely that Green’s functions on the imag¬ 
inary axis encode information in a much less usable form 
than on the real axis. 

Consider again the example of the two-site DCA for 
the single-band Hubbard model on the square lattice. In 
Ref. 39, this problem has been solved entirely on the 



FIG. 3. (Color online) Real- and imaginary-frequency Green’s 
functions computed by converging the DMFT self-consistency 
equation Eq. (11) for the two-site dynamical cluster approxi¬ 
mation to the single-band Hubbard model on the square lat¬ 
tice with next-nearest neighbor hopping t'/t = 0.3, half-band 
width D = At, interaction U = 2.5 D, and band filling n = 0.96 
in the paramagnetic phase (as in Fig. 1). See Ref. 51 for def¬ 
inition of the model and the meaning of the orbital (patch) 
quantum number a K = ±”. Panels (a-c): Electron spectral 
function A(uj) = — 4lmG ret (u;) obtained by converging on 
imaginary-frequency axis (itMPS) using number of bath sites 
and broadenings as specified in the figure, and compared to 
unbroadened (77 = 0 ) real-frequency axis computation using 
Lb/Lc = 39 bath sites per correlated site of Ref. 39. Panel 
(d): Converged Matsubara Green’s function for number of 
bath sites shown, compared to numerically exact quantum 
Monte-Carlo result of Ref. 51, computed at /3 = 200/D. 
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real axis using Lb/L c = 39 bath orbitals. Here, we con¬ 
verge the DMFT loop on the imaginary axis and compute 
the spectral function in a final real-time evolution using 
L&/L c = 3, 5, 7 bath orbitals. We compare both solutions 
in Fig. 3. Whereas for the (central) momentum patch 
“+” shown in Fig. 3(a), we find satisfactory agreement of 
the imaginary-axis with the real-axis calculation, this is 
not the case for the (outer) momentum patch ” shown 
in Fig. 3(b), even though the corresponding imaginary- 
axis Green’s function is well reproduced, see Fig. 3(d). 
Evidently, in Fig. 3(b), the central peak and the pseudo¬ 
gap at the Fermi edge are smeared out by a broaden¬ 
ing 77 = 0 . 2 D that hides finite-size effects to a large de¬ 
gree. Reducing the broadening to 77 = 0.05D as shown 
in Fig. 3(c), again reveals the pseudo-gap and the cen¬ 
tral peak; but together with unphysical finite-size effects. 
We observe that the nature of these finite-size effects is 
qualitatively comparable when using different numbers 
of bath sites Lb/L c = 3,5,7. Already for Lb/L c = 3, 
we obtain reasonable results. On the imaginary axis, by 
contrast, Lb/L c = 5, 7 still improve over Lb/L c = 3 and 
almost agree with the numerically exact QMC data for 
j3 = 200/.D of Ref. 51, see Fig. 3(d). However we empha¬ 
size that even with the modest number of bath sites used 
here the basic features of the spectral function are repro¬ 
duced (for example the areas in given frequency ranges). 


IV. THREE-BAND CALCULATIONS 
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FIG. 4. (Color online) Density per orbital as function of 
chemical potential for three-band Hubbard-Kanamori model 
Eq. (13) using the semi-elliptic density of states Eq. (14) 
and U = 12 1, obtained from single-site DMFT approximation 
evaluated using imaginary MPS (crosses) and CTQMC data 
(circles, Figure 1 of Ref. 61, inverse temperature /3 = 50 /t). 
In the DMRG computations the bath fitting was performed 
using deft = 100/t, uic = 6t and a = 1 with three bath sites 
per correlated site ( Lb/L c = 3). The maximal matrix di¬ 
mensions was m = 300 for the ground state calculation, ex¬ 
ploiting the SU(2) symmetry, which leads to the high preci¬ 
sion ((H — E ) 2 } ~ 10 -14 . For the time evolution, we com¬ 
puted G|(t) in (4a) in steps of At = 0.1 /t allowing a global 
truncation error of 5 • 10 -4 per step, up to imaginary time 
Tmai = 100/t and used linear prediction for higher times. 
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A. Three-band model in single-site DMFT 

We now demonstrate the power of the method by ap¬ 
plying it to three-band problems in the single-site ap¬ 
proximation (where comparison to existing calculations 
can be made) and the two-site approximation. Both was 
hitherto not accessible to DMFT+DMRG computations. 

We study the three-band Hubbard-Kanamori model 
with Hamiltonian (omitting the site index i in the fol¬ 
lowing definition of H\ OC i) 

H = E/ £fc C 4,a,(T C ^A<T + E H\oc,i 
k,a,b,cr i 

tfloc = ^ ^ (/^ ^a)^a,<7 

a ,<7 

+ y; Un a ^n a ^ (13) 

a 

d~ ''y ' t tl a (T Tlb^—(7 d - (U J'}t2 a a -Tlb^(j 

a>b,a 

— + dl ^dl ^d a ^d a: i + h.c.), 

a^b 

where i labels sites in a lattice and k wave vectors in the 
first Brillouin zone, ni } a tCr = d\ a a di^ a ,a is the density of 
electrons of spin a in orbital a on site i, fj, is the chemical 
potential, A a is a level shift for orbital a, e^ b is the band 


dispersion, U is the intra-orbital and U' the inter-orbital 
Coulomb interaction, and J is the coefficient of the Hund 
coupling and pair-hopping terms. We adopt the conven¬ 
tional choice of parameters, U' = U — 2 J, which follows 
from symmetry considerations for d-orbitals in free space 
and holds (at least for reasonably symmetric situations) 
for the t 2 g manifold in solids [60]. 

We study the orbital-diagonal and orbital-degenerate 
case (A. = 0) on the Bethe lattice, i.e. the non¬ 
interacting density of states is semi-elliptic, 

^,oM = ^l-(i) 2 . (14) 

In the single-site approximation, the impurity Hamilto¬ 
nian used within DMFT is given by 

H — H\ oc Hcoupl T -^bathj 
-^coupl — E vw di ,jCl, a ,a+ h.c., (15) 

l,a,a 

-^bath — ^ ^ £l,a.a < -'Z,a,cr^T c b cr ’ 
l,a,a 

where cj a a creates a fermion in the bath orbital l, V^ a ,a 
describes the coupling of the impurity to the orbital l, 
and ei t a.a denotes the potential energy of orbital l. The 
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hybridization function is then given by 


L b /L e 


E 


1 = 1 


\Vl,a, a \ 

Z - £l,a 


(16) 


Figure 4 compares the dependence of the particle den¬ 
sity n on the chemical potential /r obtained by the MPS 
methods used here to those obtained by numerically ex¬ 
act CT-QMC methods [61]. The plateaus in n(/z) are the 
Mott insulating regimes of the phase diagram. The agree¬ 
ment is very good in general, confirming the reliability of 
our new procedure even with only three bath sites per 
correlated site. This leads to an extremely cheap com¬ 
putation, for which a single iteration of the DMFT loop 
took about 30 min on two 2.8 GHz cores (see Appendix 
A 2 for more details). 

In panel (a) of Fig. 5 we show a more stringent test, 
namely the dependence of the self energy on Matsubara 
frequency, in a parameter regime where the self energy 
was previously found [62] to exhibit an anomalous w? 
frequency dependence and (in some regimes) a non-zero 
intercept as w —> 0. These phenomena are associated 
with a spin-freezing transition [61, 62], 

Panel (a) of Figure 5 shows that the known low fre¬ 
quency u> < t self energy is already accurately repro¬ 
duced even for the computationally inexpensive choice of 
Tb/L c = 3, although one observes deviations for the high- 
frequency behavior. The deviations at high frequency de¬ 
crease as the number of bath sites is increased, although 
full convergence at all frequencies has not been demon¬ 
strated. Panel (b) of Fig. 5 shows that the deviations 
are linked to the impossibility of fitting the hybridiza¬ 
tion function equally well for all frequencies using only a 
small number of bath sites. The large deviations at high 
frequencies are due to the choice a = 1 in (12), which en¬ 
forces good agreement for low frequencies and allowed to 
successfully reproduce the MIT transition in Fig. 4. In¬ 
creasing the number of bath sites to Lb/L c = 5 leads to a 
much better approximation of the hybridization function 
also for high frequencies, with concomitant improvement 
in the self-energy (Fig. 5(a)). 


B. Three-band model in two-site DCA 

We now present results obtained using a two-site DCA 
approximation to the three-band model of the previous 
section. For this problem there are no low-temperature 
results available in the literature. The size of the prob¬ 
lem is beyond the scope of standard ED. The truncated 
configuration interaction (Cl) impurity solver [20] allows 
to access a relatively high number of bath sites but is 
limited in the number of correlated sites: e.g. in Ref. 24, 
a problem with L c = 3 and Lb = 30 was computed, and 
in Ref. 56, one with L c = 4 and Lb = 20. The three- 
band two-site DCA though has L c = 6 correlated sites 
and it remains to be seen whether this is in reach for 



I_I_I_I_I_I_ 

0 0.5 1 1.5 2 2.5 3 

(On/t 


FIG. 5. (Color online) Imaginary part of Matsubara axis 
self energy E (panel (a)) and imaginary part of hybridiza¬ 
tion function A (panel (b)) for densities shown obtained from 
converged itMPS solution of single-site DMFT for three-band 
Hubbard-Kanamori model (15) for \J = 8t and J = U/6 . 
Crosses (color online) represent itMPS data and black circles 
depict CTQMC data from Figure 3 of Ref. 62, computed at 
inverse temperature /3 = 100/f. We choose all parameters as 
described in the caption of Fig. 4, in particular, for the bath 
fitting (12), we use /3 e a = 100/t, uj c = 6f and a = 1. Choosing 
a = 1 enforces agreement for low frequencies at the price of 
disagreement at high frequencies, which is observed in both 
panels (a) and (b). In panel (b), A denotes the hybridization 
function that is fitted with the hybridization function A dlscr 
of the discrete impurity model. 


the Cl solver. The problem is also challenging for stan¬ 
dard CTQMC. Recent technical improvements on mit¬ 
igating the sign problem [63] enabled Ref. 64 to treat 
this model at the temperature of T = 0.025 D with D 
the half bandwidth, which required large computational 
resources. This temperature is relatively high as the au¬ 
thors stated that in the study of a simpler two-band two- 
site model, where they reached T = 0.0125D, it was “in¬ 
tractable” to reach low enough temperatures to clarify 
whether bad metal/spin freezing behavior was intrinsic 
or not [63]. 

We study the model on the two-dimensional square lat¬ 
tice, i.e. using e^ b = —2 1 (cos k x + cos k y ) S ab . We use the 
momentum patching of Ref. 51 ; this definition is also used 
in the single-band computations of Figs. 1 and 3. We note 
that this model is not directly relevant to layered mate¬ 
rials where the t 2 g orbitals are relevant, because in the 
physical situation the two dimensionality will break the 
three-fold orbital degeneracy. However the system is well 
defined as a theoretical model and useful to demonstrate 
the power of our methods. 








FIG. 6. (Color online) Comparison of results obtained using 
MPS methods with for Lb/L c = 3 for single site (Is) and two- 
site (2s) DMFT approximations to the Hubbard-Kanamori 
model (13) on the two-dimensional square lattice with half 
bandwidth D = At, e% b = —2 1 (cosk x + cosk y ) 5 ab , U' = U — 
2 J, J = U/A and n = 3 (^ = 5U/2 — 5 J), that is, in the 
particle-hole symmetric case. Panel (a): spectral functions 
for broadening r/ = 0.2D; panel (b) broadening ?? = 0.05ZX 
In panel (c), we show the imaginary-time evolution of G > (r) 
as defined in (4a), confirming by comparison to a calculation 
for a smaller bath Lb/L c = 2 that this quantity has been 
converged with respect to the bath size. The maximal bond 
dimension for the ground state search was m = 1000. 

As is the case for the Cl method, the DMRG method 
used here is easily able to treat a large number of bath 
sites if the number of correlated sites is small: for L c = 1, 
DMRG has already often be proven to treat L f, > 120 
bath sites, and for L c = 2, L b > 80 is easily accessible 
[39, 40]. However for more correlated sites, the number of 
bath sites that can be added at given computational cost 
decreases. For L c = 6, we use L\, = 18, i.e. L&/L c = 3, 
which we have shown to be sufficient to produce reliable 
results in the single-site case without requiring overly 
large computational resources (computation time of sev¬ 
eral hours per DMFT iteration on two cores). 

We have tested the two-site calculation by converg¬ 


ing the DMFT loop for the three-band Hubbard model 
Eq. (15) with U' and J = 0 and comparing the results 
with a corresponding two-site single-band DCA. Perfect 
agreement is obtained (not shown). Non-zero values of U' 
and J create additional entanglement and make compu¬ 
tations more costly. It is then a decisive question whether 
a real-space or a momentum-space representation of the 
impurity-cluster is less entangled. We discuss this in Ap¬ 
pendix B finding that for the single-band Hubbard inter¬ 
action, both representations yield similar entanglement, 
whereas for the Hubbard-Kanamori interaction, the real- 
space representation is much less entangled. Computa¬ 
tional cost is therefore tremendously reduced by using 
the real-space representation, which comes with an off- 
diagonal hybridization function. This is the opposite be¬ 
havior as observed for QMC, where the off-diagonal hy¬ 
bridization function creates a severe sign problem. We 
further note that in the real-space representation, strong 
interactions yield a less and less entangled impurity prob¬ 
lem, as electrons become more and more localized. 

We now present results for the more physically rele¬ 
vant case U' = U — 2 J with J = U/A. For these parame¬ 
ters, at half filling the critical interaction for the MIT in 
the single-site DMFT approximation is U c ~ 1.3 D [65]. 
Fig. 6(a) shows that our results are consistent with this 
estimate: the dashed lines depict the single-site (Is) re¬ 
sults, showing a metallic solution (spectral function non¬ 
zero at w = 0) for U = D, and an insulating solution 
(spectral function zero at ui = 0) for U = 2D. In the 
two-site (2s) DCA (solid lines), by contrast, the critical 
value U c for the MIT is lowered. Even at U = D the 
ui = 0 spectral function is zero (the small non-zero value 
in panel (a) is an effect of broadening, as seen in panel 
(b)). The different nature of the metallic and insulating 
solutions is also visible on the imaginary axis in the dif¬ 
ferent nature of the decay of the imaginary-time Green’s 
function. This is plotted in Fig. 6(c) for U = D- clearly, 
a power-law decay is observed for the metallic solution 
obtained in the single-site DMFT, whereas an exponen¬ 
tial decay is obtained for the insulating solution obtained 
within the two-site DCA. 


V. CONCLUSION 

This paper introduces an imaginary-time MPS 
(itMPS) solver for DMFT and shows that it can treat 
complex models, not easily accessible with other meth¬ 
ods, at modest computational cost. This development 
establishes DMRG as a flexible low-coast impurity solver 
for realistic problems, such as encountered in the study of 
strongly-correlated materials. The crucial advance stems 
from the fact that imaginary-time evolution does not cre¬ 
ate entanglement, and hence allows to compute Green’s 
functions numerically exactly, provided a ground state 
calculation is feasible. 

The method can be improved in many ways. In par¬ 
ticular, different representations of the impurity problem 
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exhibit different degrees of entanglement, so optimizing 
the representation of the impurity problem is a promising 
route. Ideas from ED approaches for constructing rele¬ 
vant subspaces [21-24] of the Hilbert space may lead to 
further improvements. Such techniques have been suc¬ 
cessfully combined with MPS [66]. Another route to re¬ 
duce computational effort and by that reach even more 
complex models could consist in performing computa¬ 
tions for the reduced dynamics of the impurity [67], or 
making use of extremely cheap algorithms for comput¬ 
ing the Green’s function at elevated temperatures [50]. 
Finally, we note that using MPS as an impurity solver 
makes accessing entanglement as a quantity for under¬ 
standing the properties of the embedded impurity clus¬ 
ter very easily accessible. Proposals in this direction have 
been made for cellular DMFT [68] and for impurity mod¬ 
els generally [69]. 
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Appendix A: Further technical details 
1. Ground state optimization 

The main challenge in solving the ground state prob¬ 
lem of a typical cluster-bath Hamiltonian as encountered 
in DMFT, stems from the fact that DMRG is a varia¬ 
tional procedure that is initialized with a random state, 
which is then optimized locally. A local optimization 
procedure is slow when optimizing a global energy land¬ 
scape. In addition, the local optimization is prone to get¬ 
ting stuck in local minima, if no “perturbation steps that 
mix symmetry sectors” are applied. The standard per¬ 
turbation techniques for single-site DMRG [70-72] rely 
on “perturbation terms” that are produced by contract¬ 
ing the Hamiltonian with the MPS. If the Hamiltonian 
itself does not contain terms that mix the symmetry sec¬ 
tor, these methods do not work. 

A typical cluster-bath Hamiltonian has both features, 
a global variation of the potential energy and parts that 
are not connected with symmetry mixing terms, such as 
in the three-band Hubbard Kananrori model at J = 0. 
This situation is sketched in Fig. 7. 

In Ref. 39, the models under study allowed to solve 
this problem using the non-interacting solution. For the 
general models studied in the present paper, an unbiased 
numerical technique has to be employed. What we do in 


bath bath 

--cluster ^^ 

< 

-no symmetry sector mixing (e.g. density interaction) 

-symmetry sector mixing (e.g. hopping) 

E potential energy 

FIG. 7. (Color online) Sketch of a typical cluster-bath Hamil¬ 
tonian (L c = 3 , Lb = 6) when it’s mapped to a one¬ 
dimensional chain. Dashed lines depict couplings that do not 
mix symmetry sectors, and solid lines depict couplings that 
mix symmetry sectors. 

practice, is to first find the ground state of a system with 
additional symmetry-mixing couplings (denoted as red 
solid lines in Fig. 7) that are then adiabatically switched 
off. In practice, we sweep 5 to 10 times with additional 
hoppings of 10% magnitude of the physical hoppings, 
and another 5 to 10 times with additional hoppings of 
1% magnitude. After these preliminary sweeps, the quan¬ 
tum number (e.g. particle number) distribution has glob¬ 
ally converged, and we can continue with converging the 
ground state of the exact Hamiltonian. 

2. Convergence of DMFT iteration 

The calculations for the three-band single-site DMFT 
in Sec. IV A are only trivially parallelized using one core 
to compute the imaginary time evolution of each the par¬ 
ticle (>) and the hole (<) Green’s functions Ga{r). 

In Fig. 8, we show the convergence DMFT loop for the 
single-site DMFT for the three-band Hubbard-Kanamori 
model as studied in Fig. 5. Fig. 8(a) shows the con¬ 
vergence of the Matsubara Green’s function down to a 
precision of 10 -3 . Figs. 8(b) and (c) show the conver¬ 
gence of the density and of the ground state energy per 
particle, respectively. Fig. 8(d) shows the computation 
time. An iteration on the Matsubara axis takes about 
30 min. The final real-axis computation (iteration 31) is 
considerably more expensive, but can still be optimized. 

Appendix B: Least-entangled representation and 
off-diagonal hybridization functions 

1. Geometry and general considerations 

In Ref. 40, some of us showed that the star geome¬ 
try of the impurity problem can have substantially lower 
entanglement than its chain geometry. In the star geom¬ 
etry, DMRG profits from the small entanglement of the 
almost occupied states with low potential energy with 
the almost unoccupied states with high potential energy. 
A high weight for the superposition of a low- with a high- 
energy state is physically irrelevant. In the star geometry, 
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FIG. 8. (Color online) Single-site DMFT for three-band 
Hubbard-Kanamori model as studied in Fig. 5. Here for the 
case n = 1.77 (fi = 5.0) and Lb/L c = 3. To obtain the solu¬ 
tion for n = 1.79 as shown in Fig. 5, we chose fi = 5.1 and 
started from the n = 1.77 solution. Panel (a): Convergence 
of Matsubara Green’s function in the DMFT loop, starting 
from the non-interacting solution. Panel (b): Convergence 
of the density in the DMFT loop. Panel (c): Convergence 
of the ground state energy per particle in the DMFT loop. 
Panel (d): Computation time. An iteration on the Matsub¬ 
ara axis takes about 30 min. The final real-axis computation 
(iteration 31) is considerably more expensive, but can still be 
optimized. 


DMRG is able to eliminate these superpositions as po¬ 
tential energy is separated locally, i.e. in the same basis 
in which DMRG optimizes the reduced density matrix 
in order to discard irrelevant contributions. In princi¬ 
ple, as mentioned in Appendix C, ideas from basis selec¬ 
tive approaches in exact diagonalization are a different 
method to account for the fact that many states in the 
Hilbert space have a negligible weight for the computa¬ 
tion of the Green’s function and only few physically rele¬ 
vant states occupy a small fraction of the Hilbert space. 
Among these are the truncated configuration interaction 
(Cl) [20, 23, 24, 56], the basis-selective ED [21] or the 
coupled cluster methods in quantum chemistry. As these 
methods can be combined with DMRG [66], they might 
be a further route to construct efficient representations 
of the impurity-cluster problem 

In the present paper, the question of the least entan¬ 
gled representation of the impurity problem shall be re¬ 
stricted to the question of which basis to choose in a DCA 
calculation. This is of high relevance also in another con¬ 
text: In the real-space representation, the hybridization 
function becomes off-diagonal. For CTQMC, this gen¬ 
erates a sign problem. In our approach, this doesn’t af¬ 
fect computational cost much in the single-band Hubbard 
model. It even leads to a tremendous reduction of com¬ 
putational cost for the three-band Hubbard Kanamori 
interaction. 


2. DCA in momentum or real space 


The complexity of the interaction determines whether 
the real- or the momentum-space representation of the 
cluster-bath Hamiltonian is less entangled. In real space, 
the interaction has a simple form, but the hybridization 
function has off-diagonal contributions, which result in 
additional couplings of cluster and bath sites. In mo¬ 
mentum space, the hybridization function is diagonal but 
the interaction becomes off-diagonal. The additional cou¬ 
plings induced by that depend on the complexity of the 
interaction. 

Let us be more concrete. For the two-site case, the 
discrete Fourier transform yields the even and odd su¬ 
perposition of the real-space cluster. 



+ 4 ) 
^( 4 - 4 ) 


(Bl) 

(B2) 


where the index of d K labels momentum patches K and 
the index of dj labels real space cluster sites i. There 
might be further indices labeling spin or orbital. 

In real-space, the hybridization function has the form 




VjVn 


1=1 


£l ‘ 


(B3) 


where the symmetry of the real-space cluster imposes 
Aij(z) = Aji(z). In momentum space, the hybridization 
function is diagonal 


Ak(z) 


E 


i—i K —i 


V Kl V Kl 

z - ski 


(B4) 


and symmetry is reflected in the reduced number of bath 
sites per patch L' b = Lf,/i c , where L c = 2 is the number 
of momentum patches. 

We choose to use the momentum representation for 
the bath discretization, as was done for the real-axis in 
Ref. 39. While on the real-frequency axis this is the 
only viable option, the bath fitting on the imaginary- 
frequency axis via (12) is possible also for the off-diagonal 
real-space case. In real space, e.g. , particle hole symme¬ 
try can be easily imposed in the fitting procedure, while 
this is not possible in momentum space. 

Given the parameters of the momentum space repre¬ 
sentation obtained by performing a bath fit via (12), we 
define the parameters of the equivalent real-space repre¬ 
sentation as follows: In momentum space, bath param¬ 
eters are indexed by Ik = 1 L' b = Lb/L c and in 

real space, bath parameters are indexed by l = 1, ...,Lf,, 
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then 

£i = £1 ,i 1= i for l = 1, L' b 

£ i = e 2 ,i 2 =i-L' b f° r l = L' b + 1 ,..., 

V u = V 2 i = j^V Ul=l for l = 1, ...,4 

Vll = —= ^75^2 ,l 2 =l-L' b f° r ^ = -^6 + 1 ) (B 5 ) 

Whereas the momentum-space Hamiltonian has L b 
non-zero couplings Vki k , the real space Hamiltonian has 
L c x L b couplings Vu. On the other hand, the interaction 
part generates L c x (L c — 1) additional non-local couplings 
in the momentum-space representation as compared to 
the real-space Hamiltonian. 

From this one could naively expect that the real-space 
representation is less entangled if L c x (L c — 1) > L c x L b . 
Numerical experiments show that the real-space repre¬ 
sentation is much more favorable than this estimate. 
For a single-band Hubbard model, we find about the 
same entanglement in the real space and the momentum 
space representation, with slight advantages for momen¬ 
tum space. In the three-band Hubbard Kanamori model, 
the real-space representation is considerably less entan¬ 
gled and leads to a tremendous reduction of computa¬ 
tional cost. In particular, we were not able to obtain the 
results of Fig. 6 in the momentum space representation 
when using L b /L c = 3, only for L b /L c = 2 but then at 
much higher computational cost. 

Appendix C: Green’s functions from matrix product 

states 

Even though the following discussion is not needed to 
set up the imaginary-time real-time impurity solver, it 
is of general interest in this context and might stimulate 
further advancements. 

A computation of A(u) = (ipo\fi(w— {H—Eo))\ipo) via a 
computation of all eigen states of H is extremely redun¬ 
dant as only a tiny neighborhood J\f = {|-0) | (V'l ^1 V’o) 7^ 
0} of a the single-particle excitation 1 1 /> 0 ) contributes in 
the sum (inserting identities J2 n \E n )(E n \) in A{w). In 
ED, this is exploited by systematically constructing the 
subspace A f by spanning it using particle-hole excitations 
[20, 21], which might also be a viable route for further de¬ 
velopments within DMRG [66]. In DMRG, one needs to 
make a statement about the entanglement of the states 
in the subspace A f\ one might note that these are in 
general more strongly entangled than the single-particle 
excitation |^ 0 ), but should still be much less entangled 
than the rest of the Hilbertspace. This is illustrated in 
the sketch Fig. 9. 

In Ref. 30, some of us argued that expanding the spec¬ 
tral function in a family of orthogonal functions is a 
natural way to construct a basis for A/”, starting from 
the lowly entangled \i[>o) and successively increasing en¬ 
tanglement of states and thereby computational com¬ 
plexity in a sequence of basis states \ip n ). Ref. 30 dis¬ 
cussed the expansion of A(u) in Chebyshev polynomials 



FIG. 9. (Color online) Single particle excitation |gAo) of the 
ground state \Eo) and the subspace A f = {|V ) } | ^ 

0} of the Hilbert space that is relevant for the computation 
of a single-particle spectral function of the form (gAo | <5 (cu — 
H)\tpo)- The single particle excitation is very lowly entangled, 
the subspace is more strongly entangled, but still in general 
more lowly entangled than the rest of the Hilbert space. 

T n (^) = arccos(ncos(^)), which are orthogonal with re¬ 
spect to an inner product weighted by w(x) = y/1 — x 2 
[73], and in the plane waves exp (iui — ) (orthogonal with 
weight function w(x ) = 1), where the energy a is chosen 
larger than the support of A(w). The associated gener¬ 
ated sequences of basis states then are 

I iC e ) = 2(^ + 6M- 1 ) - IC-2), b G [-1,1], 

IV’n™ 6 ) = exp(—i(JJ — E 0 )^)\ipo), (Cl) 

and have different entanglement properties. The states 
l^tnne) assoc i a t e d with time evolution are in general 
less entangled than the states IV 7 )) 116 ) associated with the 
Chebyshev recursion [30]. This is due to the observa¬ 
tion that error accumulation in the Chebyshev recur¬ 
sion is worse conditioned than in time propagation [30], 
which necessitates to keep the error in a single step of the 
Chebyshev recursion much smaller than in the equivalent 
time evolution step, which in turn requires to use higher 
bond dimensions in the Chebyshev recursion making it 
less efficient. In addition to the statements of Ref. 30, 
we note here that the sequence produced by the Lanczos 
algorithm, 

|^ an > = H — a„|fc) - |V&>, Pn G R (C2) 

can be associated with an expansion of the spectral func¬ 
tion in polynomials that are orthogonal with respect to 
an inner product weighted by w(x) = A{x) [ T 4]. This is 
very efficient but numerically unstable. 

In contrast to the previous methods, which generate 
an increasingly complex basis when determining the spec¬ 
tral function to a higher and higher precision, correction- 
vector DMRG aims to optimize a state in frequency 
space, which a priori contains contributions that have 
undergone an infinitely long time evolution. As time evo¬ 
lution creates entanglement, these states are much too 
strongly entangled for an efficient treatment. They are 
“far away” from the controlled, lowly entangled single¬ 
particle excitation I'i/’o)- In order to still perform a mean¬ 
ingful computation in frequency space, one introduces 
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a so-called (Lorentzian) broadening parameter 77 that 
damps out contributions from an infinite time evolution. 
One does then not obtain the exact spectral function but 
a broadened version as in Eq. (10). The broadening pa¬ 
rameter has to be guessed a priori: If it is chosen too 
small, high entanglement prevents convergence of the cal¬ 
culation. If it is chosen too large, one will be far from the 


exact version of the spectral function. In the expansion 
methods discussed above, by contrast, one can stop the 
computation simply when it becomes too costly. If one 
has not recovered the exact A(cj) at this point, a broad¬ 
ened version can be systematically constructed with an 
a posteriori determined 77 as in Eq. (10). 


[1] W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 
(1989). 

[2] A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992). 

[3] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozen- 
berg, Rev. Mod. Phys. 68, 13 (1996). 

[4] T. Maier, M. Jarrell, T. Pruschke, and M. Hettler, Rev. 
Mod. Phys. 77, 1027 (2005). 

[5] G. Kotliar, S. Savrasov, K. Haule, V. Oudovenko, O. Par- 
collet, and C. Marianetti, Reviews of Modern Physics 78, 
865 (2006). 

[6] A. N. Rubtsov, V. V. Savkin, and A. I. Lichten¬ 
stein, Phys. Rev. B 72, 035122 (2005), arXivicond- 
mat/0411344. 

[7] P. Werner, A. Comanac, L. D. Medici, M. Troyer, 
and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006), 
arXiv:cond-mat/0512727. 

[8] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, 
M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 
( 2011 ). 

[9] M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 
(1994). 

1101 M. Capone, L. de Medici, and A. Georges, Phys. Rev. B 
76, 245116 (2007). 

[11] A. Liebsch and H. Ishida, J. Phys.: Condens. Matter 24, 
053201 (2012). 

[12] R. Bulla, T. Costi, and T. Pruschke, Rev. Mod. Phys. 
80, 395 (2008). 

[13] D. J. Garcia, K. Hallberg, and M. J. Rozenberg, Phys. 
Rev. Lett. 93, 246403 (2004). 

[14] H. Li and N.-H. Tong, (2015), arXiv:1501.07689. 

[15] P. Wang, G. Cohen, and S. Xu, Phys. Rev. B 84, 085134 
(2015), arXiv: 1410.1480. 

[16] L.-F. Arsenault, O. A. von Lilienfeld, and A. J. Millis, 
(2015), arXiv:1506.08858. 

[17] M. Schuler, C. Renk, and T. O. Wehling, Phys. Rev. B 
91, 235142 (2015), arXiv: 1503.09047. 

[18] M. Granath and H. U. R. Strand, Phys. Rev. B 86, 
115111 (2012), arXiv:1201.6160. 

[19] H. Shinaoka, M. Dolfi, M. Troyer, and P. Werner, J. 
Stat. Mech., P 2014, P0601 (2014), arXiv: 1404.1259. 

[20] G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 
186404 (2012), arXiv:1204.5783. 

[21] Y. Lu, M. Hoppner, O. Gunnarsson, and M. W. 
Haverkort, Phys. Rev. B 90, 085102 (2014). 

[22] D. Zgid, E. Gull, and G. Chan, Phys. Rev. B 86, 165128 
(2012), arXiv:1203.1914. 

[23] C. Lin and A. A. Demkov, Phys. Rev. B 88, 035123 
(2013), arXiv: 1307.4982. 

[24] C. Lin and A. A. Demkov, Phys. Rev. Lett. Ill, 217601 
(2013), arXiv:1311.5160. 

[25] K. M. Stadler, A. Weichselbaum, Z. P. Yin, J. von Delft, 
and G. Kotliar, (2015), arXiv:1503.06467. 


[26] S. R. White, Phys. Rev. Lett. 69, 2863 (1992). 

[27] U. Schollwock, Annals of Physics 326, 96 (2011), 
arXiv: 1008.3477. 

[28] U. Schollwock, Rev. Mod. Phys. 77, 259 (2005). 

[29] A. Holzner, A. Weichselbaum, I. P. McCulloch, 
U. Schollwock, and J. von Delft, Phys. Rev. B 83, 195115 
( 2011 ). 

[30] F. A. Wolf, J. A. Justiniano, I. P. McCulloch, and 
U. Schollwock, Phys. Rev. B 91, 115144 (2015), 
arXiv:1501.07216. 

[31] S. Nishimoto, F. Gebhard, and E. Jeckelmann, J. Phys.: 
Condens. Matter 16, 7063 (2004). 

[32] M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 72, 
113110 (2005). 

[33] S. Nishimoto, F. Gebhard, and E. Jeckelmann, Physica 
B: Condensed Matter 378-380, 283 (2006). 

[34] D. J. Garcia, E. Miranda, K. Hallberg, and M. J. Rozen¬ 
berg, Phys. Rev. B 75, 121102 (2007). 

[35] M. Karski, C. Raas, and G. S. Uhrig, Phys. Rev. B 77, 
075116 (2008), arXiv:0710.2272. 

[36] R. Peters, Phys. Rev. B 84, 075139 (2011). 

[37] M. Ganahl, P. Thunstrom, F. Verstraete, K. Held, and 
H. G. Evertz, Phys. Rev. B 90, 045144 (2014). 

[38] M. Ganahl, M. Aichhorn, P. Thunstrom, K. Held, H. G. 
Evertz, and F. Verstraete, (2014), arXiv:1405.6728. 

[39] F. A. Wolf, I. P. McCulloch, O. Parcollet, and 
U. Schollwock, Phys. Rev. B 90, 115124 (2014), 
arXiv: 1407.1622. 

[40] F. A. Wolf, I. P. McCulloch, and U. Schollwock, Phys. 
Rev. B 90, 235131 (2014), arXiv: 1410.3342. 

[41] C. Gramsch, K. Balzer, M. Eckstein, and M. Kollar, 
Phys. Rev. B 88, 235106 (2013), arXiv:1306.6315. 

[42] K. Balzer, F. A. Wolf, I. P. McCulloch, P. Werner, and 
M. Eckstein, (2015), arXiv:1504.02461. 

[43] K. A. Hallberg, Phys. Rev. B 52, 9827 (1995), 
arXiv:cond-mat/9503094. 

[44] T. D. Krihner and S. R. White, Phys. Rev. B 60, 335 
(1999). 

[45] E. Jeckelmann, Phys. Rev. B 66, 045114 (2002), 
arXiv:cond-mat/0203500. 

[46] S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 
076401 (2004), arXiv:cond-mat/0403310. 

[47] S. R. White and I. Affleck, Phys. Rev. B 77, 134437 
(2008). 

[48] P. E. Dargel, A. Wollert, A. Honecker, I. P. McCul¬ 
loch, U. Schollwock, and T. Pruschke, Phys. Rev. B 85, 
205119 (2012). 

[49] A. Braun and P. Schmitteckert, Phys. Rev. B 90, 165112 
(2014), arXiv:1310.2724. 

[50] D. V. Savostyanov, S. V. Dolgov, J. M. Werner, 
and I. Kuprov, Phys. Rev. B 90, 085139 (2014), 
arXiv: 1402.4516. 



13 


[51] M. Ferrero, P. S. Cornaglia, L. De Leo, O. Parcollet, 
G. Kotliar, and A. Georges, Physical Review B 80, 
064501 (2009). 

[52] M. Hochbruck and C. Lubich, SIAM Journal on Numer¬ 
ical Analysis 34, 1911 (1997). 

[53] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and 
B. P. Flannery, Numerical Recipes 3rd Edition: The Art 
of Scientific Computing, 3rd ed. (Cambridge University 
Press, New York, NY, USA, 2007). 

[54] T. Barthel, U. Schollwock, and S. R. White, Phys. Rev. 
B 79, 245101 (2009). 

[55] I. de Vega, U. Schollwock, and F. A. Wolf, (2015), 
arXiv: 1507.07468. 

[56] A. Go and A. J. Millis, Phys. Rev. Lett. 114, 016402 
(2015), arXiv: 1311.6819. 

[57] D. Senechal, Phys. Rev. B 81, 235125 (2010), 

arXiv:1005.1685. 

[58] A. Dorda, M. Nuss, W. von der Linden, and E. Arrigoni, 
Phys. Rev. B 89, 165105 (2014), arXiv: 1312.4586. 

[59] D. J. Wales and J. P. K. Doye, The Journal of Physical 
Chemistry A 101, 5111 (1997). 

[60] A. Georges, L. de’ Medici, and J. Mravlje, Annual 
Reviews of Condensed Matter Physics 4, 137 (2012), 
arXiv: 1207.3033. 

[61] P. Werner, E. Gull, and A. J. Millis, Phys. Rev. B 79, 
115119 (2009), arXiv:0812.1520. 

[62] P. Werner, E. Gull, M. Troyer, and A. J. Millis, Phys. 
Rev. Lett. 101, 166405 (2008), arXiv:0806.2621. 

[63] Y. Nomura, S. Sakai, and R. Arita, Phys. Rev. B 89,, 
195146 (2014), arXiv: 1401.7488. 

[64] Y. Nomura, S. Sakai, and R. Arita, Phys. Rev. B 91, 
235107 (2015), arXiv:1408.4402. 

[65] L. de’ Medici, J. Mravlje, and A. Georges, Phys. Rev. 
Lett. 107, 256401 (2011), arXiv:1106.0815. 

[66] Y. Ma and H. Ma, J. Chem. Phys. 138, 224105 (2013), 
arXiv: 1303.0616. 

[67] G. Cohen, E. Y. Wilner, and E. Rabani, New J. Phys. 
15, 073018 (2013), arXiv: 1304.2216. 

[68] M. Udagawa and Y. Motome, J. Stat. Mech. 2015, 
P01016 (2015), arXiv: 1406.5960. 

[69] S. S. B. Lee, J. Park, and H. S. Sim, Phys. Rev. Lett. 
114, 057203 (2015), arXiv: 1408.1757. 

[70] S. R. White, Phys. Rev. B 72, 180403 (2005). 

[71] S. V. Dolgov and D. V. Savostyanov, SIAM J. Sci. Corn- 
put. 36, A2248 (2014), arXiv:1301.6068. 

[72] C. Hubig, I. P. McCulloch, U. Schollwock, and F. A. 
Wolf, Phys. Rev. B 91, 155115 (2015), arXiv:1501.05504. 

[73] A. Weifie, G. Wellein, A. Alvermann, and H. Fehske, 
Rev. Mod. Phys. 78, 275 (2006). 

[74] W. Gautschi, Journal of Computational and Applied 
Mathematics 178, 215 (2005). 


